Nonlinear Evolution of Baryon Acoustic Oscillations 
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We study the nonlinear evolution of the baryon acoustic oscillations (BAO) in the dark matter 
power spectrum and correlation function using renormalized perturbation theory (RPT). In a pre- 
vious paper we showed that RPT successfully predicts the damping of acoustic oscillations; here we 
extend our calculation to the enhancement of power due to mode-coupling. We show that mode- 
coupling generates additional oscillations that are out of phase with those in the linear spectrum, 
leading to shifts in the scales of oscillation nodes defined with respect to a smooth spectrum. When 
Fourier transformed, these out of phase oscillations induce percent-level shifts in the acoustic peak 
of the two-point correlation function. We present predictions for these shifts as a function of red- 
shift; these should be considered as a robust lower limit to the more realistic case that includes in 
addition redshift distortions and galaxy bias. We show that these nonlinear effects occur at very 
large scales, leading to a breakdown of linear theory at scales much larger than commonly thought. 
We discuss why virialized halo profiles are not responsible for these effects, which can be understood 
from basic physics of gravitational instability. Our results are in excellent agreement with numeri- 
cal simulations, and can be used as a starting point for modeling BAO in future observations. To 
meet this end, we suggest a simple physically motivated model to correct for the shifts caused by 
mode-coupling. 



I. INTRODUCTION 

The imprint of baryon acoustic oscillations (BAO) in 
the power spectrum and two-point correlation function of 
galaxies has emerged as a potentially robust probe of the 
cause of the present acceleration of the universe, through 
their sensitivity to the angular diameter distance and the 
Hubble constant as a function of redshift. Measurements 
of this signature can be used to help probe dark energy or 
large-distance modifications of gravity as the explanation 
for the acceleration of the universe. 

The basic idea behind the BAO method is to use the 
characteristic size of the sound horizon imprinted in the 
spatial correlation properties of galaxies or other trac- 
ers as a standard ruler In the power spectrum, 
this characteristic scale and its harmonics lead to os- 
cillatory features, the so-called "baryonic wiggles" 
In the two-point correlation function, these features 
translate to a broad peak at the sound horizon scale, 
about 100/i _1 Mpc for acceptable cosmological param- 
eters. The BAO have been detected recently in both 
the correlation function and power spectrum of galax- 
ies @,g. 

This acoustic signature evolves with time, and the key 
issue is how. According to linear perturbation theory 
the acoustic signature increases in amplitude but its spa- 
tial pattern remains static, i.e. the characteristic scale 
imprinted in the early universe remains unaltered. How- 
ever, because gravitational instability is a nonlinear pro- 
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cess, the same motions that lead to the growth of corre- 
lations also change their shape. Thus one must check the 
validity of linear perturbation theory arguments, even at 
scales as large as 100 hr 1 Mpc. This is even more so given 
the precision attainable in forthcoming experiments, e.g. 
a shift in the acoustic scale of one percent generates sys- 
tematics in the deduced dark energy equation of state 
parameter w of about five percent [3, |5|, which is not 
negligible compared to the expected statistical errors in 
the next generation of galaxy surveys. 

Previous work on the subject has been based on 
numerical simulations d, H, 0, H, Q, using the halo 
model flo| . or combinations of analytic and numerical 
techniques [H M, GS El US 03] ■ The picture emerging 
from this body of work regarding systematic effects due 
to nonlinear clustering is still far from converging into 
a coherent framework. For example, [Io| show from the 
halo model that any clustering-induced systematics must 
be small, and [l3[ argue that any systematic shift (i.e. not 
related to random motions) of the acoustic scale for dark 
matter in real space must be of unobservably small, of or- 
der 0.01% at z = 0. On the other hand, [ijj showed that 
the nonlinear power spectrum fitting formula of [l7| pre- 
dicts percent-level shifts in the acoustic peak of the two- 
point function. Similarly, [l4| argue from studying the 
power spectrum of halos in simulations and perturbation 
theory that significant shifts are expected. In addition, 
Q do a detailed study of the impact of nonlinearities on 
the power spectrum of galaxies from semianalytic models 
to conclude percent-level shifts cannot be excluded. Fi- 
nally, [13] presents evidence from numerical simulations 
that the acoustic peak of dark matter and their halos do 
experience non negligible systematic shifts, which cannot 
be explained by random motions alone, and construct a 
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physical model based on the coherent infall of pairs to 
understand their origin. 

This state of affairs is perhaps not too surprising given 
that the effects involved are small, and require great accu- 
racy from analytic and numerical methods. In this paper 
we consider this issue by using renormalized perturba- 
tion theory (RPT [l8|, [Hj]), a new approach to follow 
nonlinear clustering that includes in a systematic way 
all nonlinear effects in the fluid approximation around 
a given scale (52j. Here we concentrate on fundamental 
questions such as 1) can nonlinear effects generate shifts 
in indicators of the acoustic scale large enough that may 
bias determinations of cosmological parameters? 2) if 
so, what physics is responsible for this? Is it related 
to large-scale nonlinearities that we can hope to model 
accurately, or more complicated physics related to virial- 
ized dark matter halos? We shall see that the answer to 
the first question is 'yes', and the answer to the second 
question involves large-scale physics, that we discuss in 
detail. Our discussion emphasizes the shifts generated 
by mode cou plin g, which constitutes a new result (see 
also pjj). In |l9l ] we have already discussed in detail the 
effects of random motions in terms of large-scale physics; 
we briefly discuss these here as well in more accessible 
terms. That large-scale random motions are responsible 
for the damping of the linear power spectrum has also 
been recognized in d, [TH, |22j . 

In the present paper we concentrate on predictions 
from RPT for the power spectrum and the two-point 
correlation function. A detailed account of the tech- 
nicalities involved in calculating two-point statistics in 
RPT and their comparison with numerical simulations is 
left for a separate publication [23|. Here we present the 
main results regarding BAO for dark matter in real space 
and discuss how RPT can shed some light on practical 
parametrizations of these nonlinear effects in a more gen- 
eral situation when redshift distortions and galaxy bias 
are also present. No familiarity with RPT is assumed, the 
main ideas behind RPT and results on two-point statis- 
tics are explained in simple terms in the following section, 
while the analytic expressions for the power spectrum are 
presented in Appendix [Al 



II. RPT AND TWO-POINT STATISTICS 



Basics of RPT 



Standard perturbation theory (PT, see [24| for a re- 
view) is an expansion of the equations of motion around 
their linear solution, assuming fluctuations are small. 
Schematically, for the power spectrum this expansion 
reads 



P{k,z) = D 2 + {z)P (k) + (k, z) + P21oop(fc, z) 



where D + (z) is the growth factor at redshift z, Po(k) 
is the initial power spectrum (at high redshift) so that 
linear evolution reads Pu n (k,z) — [D + {z)} 2 P^{k) . In 
Eq. ©, 

-Plloop ^ OCfltoAun), P 2 loop ~ 0(fl in Af in ), 
and so on, where Ai; n = 47rfc 3 Pij n measures the ampli- 
tude of fluctuations at scale k in linear theory. For scales 
approaching the nonlinear regime where Au n > 1, trun- 
cation at any finite order in PT is not meaningful, as 
neglected higher-order contributions are important. 

In Renormalized Perturbation Theory (RPT, [HI), the 
main idea is to get around this limitation of PT, by mak- 
ing a resummation of an infinite subset of contributions 
to the PT expansion. As a result of this process of resum- 
mation, where terms of different order have been grouped 
together into physical objects, what remains is a new se- 
ries expansion which is not a perturbative expansion in 
the amplitude of fluctuations and, most importantly, ex- 
hibits a very different behavior: truncation at finite order 
in RPT does take into account all nonlinearities from the 
largest scales down to a given scale, the impact of smaller 
scales described by the neglected terms is highly sup- 
pressed. One the main insights that follows from RPT is 
that if we write the growth factor as, 



D+(z) = 



(4n(k,z)Mk')) 

(5o(k)5o(k')) ' 



(2) 



where S denotes the density contrast, and 6n n (k, z) — 
D+(z) <5o(k) is linear evolution (with D+ = 1 at the ini- 
tial condition); then a whole set of nonlinear contribu- 
tions to Eq. |T]) (or any correlation function) effectively 
"renormalize" the growth factor to the following, fully 
nonlinear quantity 



D + {z) — ► G(k,z) 



(5(k,z)S (k')) 
(5 (k)S Q (k')) ■ 



(3) 



where S(k, z) is the fully nonlinear density contrast. The 
function G(k,z) is known as the propagator, which can 
be thought as a measure of the memory of initial condi- 
tions, since it gives the time "propagation" of the cross- 
correlation between initial and final density contrasts, 
(S(k,z) 5 (k')) = G(k,z)(S a {k)S {k')). Note that this 
property means that all the terms in Eq. ([T]) that are 
proportional to Po (including those in the loop contri- 
butions) are resummed into G 2 Po, whereas in the re- 
maining loop terms the time dependence is dictated by 
the propagator instead of the growth factor, which essen- 
tially means to use Eq. ([3]) to replace the linear propaga- 
tion in between nonlinear interactions that make up the 
loop contributions [531 ] . 

The asymptotics of the propagator are easy to under- 
stand: at large scales, linear perturbation theory becomes 
a good approximation and thus 



G(0,z) = D + (z). 



(4) 



(1) 



On the other hand, at small scales where nonlinear ef- 
fects are dominant the cross-correlation must be driven 
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to zero, as the final density field resembles very little 
what it was at the beginning. Thus, we expect on phys- 
ical grounds that 



G(k, 







(5) 



It is this last property that is impossible to capture at a 
fixed order in PT, where G becomes arbitrarily large in 
this limit (positive or negative, depending on truncation 
order], thus leading to unphysical results. As we showed 
in [191 ] , in the high-/c limit the dominant behavior of the 
propagator can be calculated exactly giving 



G(k,z) a D+{z) exp 



D+(z) 



1 



(6) 



for ka v D + 3> 1, where the characteristic scale of decay 
determines the breakdown of linear perturbation theory 
and it is given by 



1 



r 



(7) 



i.e., the effective one-dimensional amplitude of large-scale 
velocity flows. Due to the shape of the CDM spectrum 
of fluctuations, this scale is rather large, at z = 



[D+(0)a v ] 



0.15 hMpc~ 



(8) 



while e.g. at z = 1,3 this becomes 0.24, 0.46 /iMpc , 
respectively. When kcr v D + < 1 the expression for the 
propagator is much more complicated, but it is still 
well approximated by a Gaussian with a weakly scale- 
dependent width. See [l9[ for a comparison of the RPT 
propagator against numerical simulations. 

These results have two immediate consequences [15, 
Il9j |. First, the validity of linear perturbation theory is 
much more restricted than commonly thought, as we dis- 
cuss in more detail below. Second, when such a renor- 
malization is performed, the resulting perturbation the- 
ory has the property that large scales are effectively 
"shielded" from small scales, e.g. Eq. © suggests that 
the highly nonlinear regime where the fluid approxima- 
tion breaks down has an exponentially small impact on 
large scales, since the scale in Eq. © is so large com- 
pared to scales where more complicated physics is likely 
to enter. This makes RPT a robust method for probing 
nonlinear effects at relatively large scales such as those 
involved in BAO. 



B. Two-Point Statistics in RPT 

When contributions to Eq. (JTJ are absorbed into the 
renormalization of the growth factor as in Eq. ([3]), the 
equation for the power spectrum in RPT reduces to an 
exact but simple form, 



P(k,z) = G 2 {k,z)P {k)+P MC {k,z) 



(9) 



where all contributions which are proportional to the ini- 
tial spectrum of fluctuations are now included in the first 



term, which represents how much of the primordial power 
remains at a given scale after nonlinear evolution, and 
thus has direct information on the linear power spectrum. 
Note that G describes how a given k mode propagates 
in time in the presence of all the other Fourier modes, 
and thus it depends on the initial power spectrum at 
all scales, as it is manifest in Eqs. (JSC])- By contrast, 
the second term Pmc represents the power generated by 
mode-coupling at smaller scales, and depends on the lin- 
ear power at different scales than k through complicated 
convolutions. This has a similar loop expansion to that 
sketched in Eq. (JTJ) , 



Pmc(K z) = pU°°*(k, z) + P*°° p (k, z) + . 



(10) 



but with the important difference that its time depen- 
dence is also dictated by the propagator. As a result 
of this, the convergence properties of the remaining per- 
turbation theory for Pmc are drastically changed: only 
a few terms are needed to describe nonlinear effects at 
a given scale (l8| . This is clear from the physical in- 
terpretation of each term in Eq. (fT0"|) . A term with n 
loops describes the effect of {n + 1) modes coupling to 
the mode of interest, i.e. k = ki + . . . + k„ + i, and thus 
Piioop ~ 0(A n ) while in addition there is an exponen- 
tial suppression at high-fc because the time dependence 
of P n ioop is built out of G's. Thus each successive term 
is increasingly suppressed at low-fc where A <C 1, then 
peaks at a higher k where A > 1, before decaying when 
G — ► 0. See Section UlI Bl for more discussion on this, and 
also Fig. [5] below for how the first two terms in Eq. (fTT 
describe the large-scale properties of Pmc ■ 
The two-point function follows from Eq. ©, 



f(r,z) = [G 2 ®Zo](r,z) + £ M c(r,z), 



(11) 



where the symbol <8> indicates a convolution. Since in 
Fourier space the propagator G is approximately Gaus- 
sian, the first term convolves the initial correlation func- 
tion with an approximately Gaussian kernel, leading to 
a smoothing of any features present in £o> that also in- 
duces shift of the peak location due to the fact that the 
linear correlation function is not quite symmetric about 
the peak. The second term, due to mode-coupling, natu- 
rally leads to a shift of the acoustic peak in the two-point 
function, as we will discuss in detail in Section IV Bl 

Figure Q] illustrates these results at z = 0, for the 
power spectrum (left panel) and the two-point correla- 
tion function (right panel). The acoustic signature in 
the linear spectrum consists of a roughly harmonic se- 
quence of peaks and troughs exponentially damped due 
to Silk diffusion. The damped sequence extends up to 
k ~ 0.2 ft Mpc -1 while the first trough is around k « 
0.05ftMpc _1 . In the correlation function the acoustic 
signature is a single broad peak at the sound horizon scale 
(shown also as a vertical dashed line, at 106 ft -1 Mpc for 
the cosmology we consider here) , and its width is related 
to the decreasing envelope of oscillations in the power 
spectrum 
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k [h/Mpc] r [Mpc/h] 

FIG. 1: The power spectrum (left) and two-point correlation function (right) as described by RPT, following Eqs. © and (jlll) . 
respectively. Each contribution is labeled in the same line styles left and right, and these are just Fourier transforms of 
each other. The damping of the linear spectrum by G 2 (left) translates into a smoothing of the acoustic peak (right), whose 
unperturbed position is at 106 Mpc in the cosmological model we use in this paper (dashed vertical line). The mode 
coupling power shows a strong scale dependence (left) and oscillations (of amplitude too small to be seen in this figure, see 
Fig. [6]l and this translates into a shift towards small scales in correlation function space (right). Note that £mc is positive left 
of the peak (solid) and negative right of the peak (dashed) . 



Figure Q] clearly shows the ingredients of RPT at work, 
following the decomposition in Eqs. (0 and (fTTj) . The 
left panel shows the damping of the linear features from 
G 2 Po as nonlinear motions decorrelate the Fourier modes 
with respect to their initial values, and the sharp rise 
of the mode-coupling contribution Pmc describing the 
newly generated power that dominates at small scales. 
One important conclusion from this is that the validity 
of linear theory is much more restricted than commonly 
thought [HI, OH : the comparison between linear and non- 
linear power gives a misleading picture. Deviations from 
linear theory of order 10% in power happen at scales as 
large as k — 0.05 ft. Mpc -1 , and e.g. at k = 0.13 h Mpc -1 
where linear and nonlinear power differ by less than 2%, 
RPT tells us that approximately only half of the power 
at that scale is related to the primordial power at the 
same scale, the other half is coming from other scales 
(see Fig. [3] below for more discussion on this). Note in 
addition that this compensation in power between the 
decay of G and the growth of Pmc is only approximate 
since these terms depend differently on the initial power 
spectrum and thus cosmological parameters. 

In real space (right panel), the decomposition of non- 
linear clustering effects in Eq. (|11[) is very helpful in pro- 
viding a qualitative understanding of the cosmological 
evolution of the acoustic signature. The smoothing of 



the peak by convolution with G 2 can be understood from 
the Fourier picture by the fact that a peak in real space 
requires "well prepared" Fourier coefficients, and as these 
are decorrelated by non linear evolution the peak de- 
cays in amplitude compared to linear amplification. This 
decorrelation can be understood in real space as resulting 
from the transport of matter from one point in space to 
another by large-scale flows: indeed, the result in Eq. ([5]) 
was obtained by resumming exactly the effects of the v- V 
terms in the equations of motion [19j . This explains why 
the characteristic scale of decay in the propagator is the 
variance of the velocity field, Eq. ([7]). 

Finally, note that £mc changes sign at the acoustic 
scale, leading to a shift in the peak of the two-point func- 
tion to smaller scales. This is due to the fact that Pmc 
contains information about the acoustic scale. Indeed, as 
we show in section [Vj Pmc contains acoustic oscillations, 
although of an amplitude too small to be visible in Fig. [1] 



III. RPT VS. NUMERICAL SIMULATIONS 
A. Numerical Simulations 

In order to quantify the accuracy at which RPT can 
describe the evolution of BAO we now proceed to test it 
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against measurements in N-body simulations. A precise 
description of BAO is somewhat challenging for numeri- 
cal simulations because this signature lies at large scales 
of order ~ 100 ft. -1 Mpc. Therefore a very large simu- 
lation volume is required to reduce the cosmic variance 
to levels where it can be dissentagled from true nonlinear 
evolution, also allowing the oscillations in the power spec- 
trum or the peak in the correlation function to be well 
resolved. Hence we run a set of 50 N-body realizations in 
a cubic volume of side L^ ox — 1280 h~ x Mpc and N par — 
640 3 particles with outputs at z = 1, 0.5, 0. This gives us 
a total sampled comoving volume of ~ 105( h^ 1 Gpc) 3 . 
Initially in this project we used a smaller set of 8 real- 
izations with volume (1024 h~ x Mpc) 3 and N par = 512 3 
particles with outputs at z = 2,1,0.3,0, which we only 
use here in Fig. [5] to cover a larger range in redshift. 

The initial conditions were set at z% = 49 using 
the public 2nd order Lagrangian Perturbation Theory 
(2LPT) code described in [25J, |26j, and the subsequent 
evolution was followed using the Gadget2 code [271 ]. Note 
that we used the transfer function output at z — 49 in- 
stead of z — 0, this means our acoustic peak is somewhat 
smaller in amplitude than it should be. We comment 
on this issue in section IV Bl The cosmological parame- 
ters were set to TO = 0.27, £!a = 0.73, Sit = 0.046 and 
h = 0.72. The initial power spectrum had a scalar spec- 
tral index n s = 1 and was normalized to give as = 0.9 
when linearly evolved to z = 0. We measured the prop- 
agator, power spectrum and correlation function at var- 
ious redshift outputs. In all cases the errors reported in 
this paper are for the mean of the ensemble, i.e. corre- 
sponding to a volume of ~ 105( h" 1 Gpc) 3 , and obtained 
from the scatter among the 50 realizations. Except in 
Figure [2] where we use the smaller set of simulations, with 
8 realizations giving a total volume of ~ 8( h~ x Gpc) 3 . 

In this paper we focus on a single cosmological model. 
A more detailed assessment of RPT and numerical simu- 
lations, including different cosmologies, will be presented 
elsewhere . Note however than in [l9[ we validated 
the RPT predictions for the propagator for different red- 
shifts and both densities and velocities. Since the evolu- 
tion of CDM perturbations is not self-similar, this was 
a non-trivial test. Here we concentrate on the den- 
sity power spectrum at different redshifts. Higher-order 
statistics will be considered in a forthcoming paper (29l | . 



B. The Power Spectrum 

We now present results for the evolution of the baryon 
wiggles in the power spectrum using RPT and compare 
them to numerical simulations. According to Eq. §§§ 
one needs to calculate the nonlinear propagator G and 
the mode coupling power Pmg- Here we explain the ba- 
sic results in words and refer the reader to Appendix [5] 
for more details about the calculation; a throughout dis- 
cussion of the technical aspects involved will presented 

in m. 



In [H[ we obtained an analytical prescription for the 
nonlinear propagator. This was done as follows. First, 
we calculated its \ow-k behavior from one-loop PT and 
its large-fc limit, Eq. ©, by resumming the infinite sub- 
set of diagrams that provided the dominant contribution. 
We then noticed that these two limits can be matched in 
a unique way, and without introducing any free parame- 
ters, if one regards the low-fc limit as an expansion of the 
Gaussian in Eq. (J6j) . This procedure resulted in a prop- 
agator, for both density and velocity divergence fields, 
that was shown to be in very good agreement with sim- 
ulations for all times and scales. 

The mode coupling power is decomposed as an infi- 
nite sum of partial contributions, Eq. (|10p. each of which 
is dominant in a narrow range of scales and negligible 
otherwise. Roughly speaking, as discussed above, each 
contribution is a positive "bump" centered at increas- 
ingly higher k [l8[ , see Fig. [5] below for the first two con- 
tributions. Since the BAO do not span a large range 
of scales, only few terms in the mode-coupling power 
are needed to describe them accurately. Indeed, as we 
shall see in more detail below, the first two contribu- 
tions describing two and three-mode coupling suffice; i.e. 

P M C — -Pmc° P + ^mc OP is au we nee d m Eq. ((TUJl . 

Figure [5] shows our results for the dark matter power 
spectrum as a function of scale for redshifts z = 0, 0.3, 1, 2 
in four panels (as labeled). The symbols with error bars 
correspond to measurements in 8 realizations of a cubic 
box of side 1024 h~ x Mpc. In addition we also show the 
predictions from standard one loop PT [30] | and from the 
nonlinear power spectrum fitting formula halof it [l7l |. 
The dashed line shows the linear power spectrum. Note 
that all quantities have been divided by a smooth BBKS 
power spectrum [28| with shape parameter T — 0.15 to 
increase the contrast of the acoustic signature. 

We can see that the RPT prediction matches the mea- 
surements from the simulations very well at all redshifts. 
Halofit underestimates the power at all scales and be- 
comes a better fit at higher redshift where nonlinearities 
are weaker. At z = the underestimate is about 4% at 
k = 0.1/iMpc -1 and 8% at k = 0.2/iMpc _1 . As is well 
known, the prediction from one-loop PT overestimates 
the power spectrum, at z = by about 4% (18%) at 
k = 0.1 /iMpc -1 (0.2/iMpc _1 ). Although this prediction 
improves at higher redshifts it never becomes accurate, 
e.g at z = 1 and k = 0.2 h Mpc -1 the overestimate is of 
order 7%. This statement might be in conflict with the 
conclusions in [ll[, where a better than 1% agreement is 
claimed at these scales and z > 1. The difference may 
be due to our higher as value that makes nonlinearities 
stronger. 

The measurements presented in Fig. [5] show an 
anomaly at k ~ 0.07 h Mpc -1 , where the measured power 
is significantly below the predictions. We can trace 
this to the initial conditions of one of the eight real- 
izations. This shows that even in such a large volume 
~ 8( Gpc) 3 the mean power may still differ from the 
expected cosmic mean at BAO scales. In what follows 
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FIG. 2: Nonlinear evolution of the acoustic oscillations in the dark matter power spectrum. In all cases we show the nonlinear 
power spectrum divided by a smooth spectrum 28] to make the acoustic oscillations more visible. The square symbols with 
error bars correspond to measurements in N-body simulations, whereas RPT prediction is represented by a solid red line as 
labeled. One-loop Perturbation theory (solid black line), halof it (solid magenta line) and linear theory (dashed blue line) are 
also shown. The different panels correspond to z — 0, 0.3, 1, 2 (top left, top right, bottom left and bottom right respectively). 
The agreement between the RPT prediction and the N-body measurements is excellent for all redshifts, see Fig. [3] for a more 
detailed comparison. 



we use 50 realizations of 640 3 particles in a cubical box 
of side 1280/i^ 1 Mpc that cover a much larger volume 
~ 105( h" 1 Gpc) 3 and effectively eliminates this issue. 
Based on the smaller set of simulations 14] claim that 
the suppression of the nonlinear power spectrum with re- 
spect to linear theory is of order 5% at z = 0. Based on 
the much larger volume simulations, we see a maximum 
suppression of 3.5% at k = 0.07 ft Mpc -1 for the same 



cosmology. 

Figure [3] shows a different normalization than Fig. [2j 
this time dividing by the nonlinearly evolved smooth 
power spectrum, thus taking out the scale dependence in- 
duced by the mode-coupling power. This allows us to re- 
duce the vertical scale of the plots and appreciate in more 
detail the comparison between RPT and simulations. We 
show linear theory (short-dashed line) and the different 
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FIG. 3: As Fig[2] but here we divide by the nonlinear counterpart of P smoo th to take out the scale dependence induced by 
mode coupling and thus provide a more detailed comparison of RPT predictions against simulations. Short-dashed lines denote 
linear theory, long-dashed lines denote G 2 P , thin solid lines G 2 P + -Pmc° P > and thick solid lines G2p o + p mc° P + p mc° P - For 
k > 0.2/tMpc _1 at z = and k > 0.3/iMpc _1 at z = 1, the missing higher than two-loop contributions in the mode-coupling 
power, Eq. (|10p . become important. RPT and simulations agree to better than 1%. 



contributions that enter into our RPT calculation: G 2 Pq 
(long-dashed line), G 2 P + P^ op (thin solid line), and 

our most complete calculation P ~ G 2 Pq+P^? p +P^? P 
(thick solid line). For k > 0.2 Mpc" 1 at z = 
and k > 0.3 /i Mpc -1 at z = 1, the missing higher 
than two-loop contributions in the mode-coupling power, 
Eq. (jlpp . become important. The nonlinear smooth 
power was computed at the same level of approxima- 
tion, P = G 2 P smooth + Pmc° P + P MC P > where for the 
propagator and mode-coupling power we used those of 
the initial spectrum with wiggles. This minimizes the 
impact of small numerical integration errors that may 
affect slightly differently the smooth and baryonic spec- 
trum. In addition, these ingredients depend weakly on 
detailed features as they involve integrations over broad 
ranges of wavenumbers. In detail, though, nonlinear cor- 
rections are sensitive to features in the linear spectrum, 
see discussion in section fVl where we show that the mode- 
coupling power contains oscillations. We postpone until 
then a discussion of the shifts in the acoustic pattern of 
the power spectrum. 

From this comparison we see that our calculations 
agree with simulations to better than 1%. At such level 
of detail, however, some of our approximations within the 
framework of RPT are likely to matter, see Appendix [A"l 
Similarly, as far as we know simulations have not been 
tested to this level so far (see e.g. [3ljj ) . Note that test- 



ing at this level at these scales requires large simulation 
boxes, for example according to RPT at k = 0.2 h Mpc" 1 
a simulation in a box of size 500 h~ x Mpc will have a sys- 
tematic error of 1% due to its finite volume. The level 
of agreement in Fig. [3] is nevertheless very encouraging. 
Note in particular that the thin-solid line (G 2 Po + P^°? p ) 
is a good approximation to the behavior of the first oscil- 
lation, whereas the addition of the two-loop contribution 
helps to describe the higher harmonics. This is important 
for the description of the acoustic peak in the two-point 
function, which is mostly sensitive to the fundamental 
oscillation. 



From Figs. [2] and [3] the two effects described in Eq. ([9]) 
are clearly visible. At lower redshifts more nonlinear 
evolution increases the damping of the acoustic signa- 
ture and also increases the scale-dependence induced by 
mode-coupling. These deviations from linear theory are 
substantial, as discussed in Section IIIBI above. For ex- 
ample, at z = the correction to linear evolution at 
k = 0.05 h Mpc -1 is already about 10% in power, and 
while at k = 0.13 h Mpc" 1 linear and nonlinear power 
agree to better than 2%, only about half of the nonlinear 
power is coming from power present in the initial spec- 
trum at the same scale, i.e. G 2 Po ~ P/2. 
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FIG. 4: The two-point correlation function for z = (left panel) and z = 0.5 (right panel). The dashed line corresponds to 
the linear two-point function, the solid line is the prediction of RPT and the symbols with error bars are the measurements in 
numerical simulations, corresponding to 50 realizations comprising a total volume of 105 (ft -1 Gpc) 3 approximately. 



C. The Two-Point Correlation Function 

To calculate the two-point correlation function in RPT 
we Fourier transform the RPT power spectrum predic- 
tion presented in the previous subsection, although here 
we only use the one-loop approximation to the mode- 
coupling power, i.e. we transform P ~ G 2 Pq + P^c 9 
to real space. The two-loop contribution to Pmc is not 
included since it only introduces very small corrections 
at BAO scales, and performing its Fourier transform re- 
quires a very accurate evaluation, which is numerically 
costly. 

Figure 2] shows the prediction of RPT for the two-point 
correlation function (solid line) against the measurements 
in N-body simulations (symbols with error bars) and 
the linear theory correlation function (dashed line) for 
a broad range of scales. The left panel shows z = and 
the right panel corresponds to z = 0.5. 

The agreement between RPT and N-body measure- 
ments for the two-point function is remarkable, although 
expected from the results on the power spectrum pre- 
sented in Figs. [2] and [3] The different actions of G 2 and 
£mc m real space is another way of seeing that the can- 
cellation in the power spectrum between G 2 Pq and Pmc 
presents a somewhat misleading picture. The action of 
these two effects is completely different in the correlation 
function, and one clearly sees large (~ 30%) deviations 
from linear theory at 100 ft. -1 Mpc scales. This is because 
measuring the power spectrum at a given scale doesn't 



say how much of it is correlated with the initial con- 
ditions and how much is due to mode-coupling, unless 
one also measures G, something one cannot do in ob- 
servations but it is simple enough to do in simulations 
(see [HHil and Section [TV] below) . As we demonstrate 
in section IV C[ the mode-coupling power leads in correla- 
tion function space to contributions which involve deriva- 
tives of the linear correlation function, and when features 
are present these terms can become important at large 
scales. Similarly, convolution with G 2 becomes a signifi- 
cant nonlinear effect when the linear correlation function 
has features which have a width comparable or smaller 
than that of G 2 (given by 2a v ), which from Eq. ([5]) is 
~ 13, 8, 4 h^ 1 Mpc at z = 0, 1, 3. The acoustic signature 
satisfies this condition at low rcdshift. 

We postpone discussion of the evolution of the peak of 
the two-point function until section IV Bl where we study 
the shift of the acoustic peak as a function of redshift. 



IV. RPT AND THE HALO MODEL 

Figure [TJ and Eqs. © and (TTJ shows that RPT 
presents a very similar decomposition to that in the halo 
model (see [32[ for a review). In this case the density 
field is modeled as a collection of halos containing all the 
mass 

p(x) = ) mj u mi (x - Xj), (12) 



0.01 



0.05 0.1 
k [h/Mpc] 



FIG. 5: The mode-coupling power spectrum Pmc, see Eq. (|20|l . measured from N-body simulations (symbols with error bars) as 
a function of scale at z — 0. The horizontal line gives the prediction of the halo model, the 1-halo term asymptotics in Eq. QB) , 
valid at the scales shown in this plot. The other lines give the prediction of Pmc in RPT in one and two- loop approximation 
(as labeled) and their sum (top solid line). For k > 0.2/iMpc -1 higher than two-loop contributions not included here become 
important. Note that the low-fc asymptote k 4 is achieved as predicted by momentum conservation. The halo model grossly 
violates this behavior. Although not apparent in this figure, the mode-coupling power contains oscillations, see Fig. [6] 



where u m denote the profiles of halos of mass m, i.e. 
how the mass inside each halo is spatially distributed. 
The density power spectrum can written as 

P{k) = P lh (k) + P 2h (fc) (13) 

where the "one-halo" term represents the contribution 
of objects that are in the same halo (thus dominating 
at small scales), and the "two-halo" term dominates at 
large scales and represents the contribution of pairs in 
different halos. They are given by, respectively 

P lh (fc)=(2^) 3 J dm(j\ n(m)\u m (k)\ 2 , (14) 

P 2h (k)=[b 1 (k)) 2 P ha (k), (15) 
with a scale-dependent bias factor 

bi(k) = (2tt) 3 J dm f^j \ n(m) bi(m) u m (k), (16) 



where n(m) is the mass function, &i(m) is the linear bias 
factor for halos of mass m and u m (k) is the Fourier trans- 
form of the halo profile, normalized given Eq. (|12|) so that 

(2tt) 3 u m (k = 0) = J d 3 x u m (x) = 1. (17) 

In [l8| we showed that interesting similarities can be 
established between the halo model and RPT, in par- 
ticular, the scale dependent bias factor in Eq. (fTS"]) was 
shown to be b\{k) = Gh(k, z)j D+(z) where Gh is the 
propagator calculated in the halo model following the 
rules of RPT. Written in this simplified and often used 
form the two-halo term would correspond to propagator 
renormalization in RPT, i.e. the term G 2 Pq in Eq. ([9]). 
The one-halo term can then be identified with the mode 
coupling power since it describes mode-coupling in the 
same sense as Pmc i n RPT: these terms represent the 
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contribution to the power that is not proportional to the 
linear spectrum at the same scale. 

Here we would like to point out some important differ- 
ences between RPT and the halo model as written above, 
which point to serious shortcomings of the latter for BAO 
modeling. First, in the halo model as written above the 
characteristic scale of decay of BAO is related to the virial 
radius of halos through the dependence on the halo profile 
in Eq. (|16p . and the resulting spatial scale is much smaller 
than that from Eq. ([5]) thus "delaying" the decay until 
higher k's. Furthermore, the decay as a function of scale 
is only power-law, instead of Gaussian. All this leads to 
a halo-model propagator or bi(k) that decays too slowly 
as fc increases. This is not too surprising, as RPT tells 
us that large-scale nonlinear effects unrelated to smaller 
scale virialization physics are responsible for the decay of 
the propagator. The obvious way to improve this situa- 
tion would be to add nonlinear contributions from both 
dynamics [IS [34], HBj] and bias [3, to the two-halo 
term in order to try to model BAO scales. 

The second, related, and more difficult to fix problem 
is that the one-halo term has a very different behavior 
at large scales than the mode-coupling power in RPT. In 
fact, in the low-fc limit appropriate for BAO scales using 
Eq. |T71) in Eq. (O it follows that 



Pih(0) = 



dm 



— n(m) = 1.96(/i _1 Mpc) 3 
P 



where we have used the Sheth-Tormen mass function [371 ] 
at z = 0. This means that on scales relevant to BAO the 
halo model power spectrum, Eq. (|13p . reads, 



P(k) ~ [h(k)] 2 P lin (k) + c, 



(19) 



where c = Pih(0). This leads to a very different pic- 
ture of the impact of nonlinearities on the acoustic sig- 
nature. Indeed, Eq. (|T9|) suggests that while in power 
spectrum space the baryonic wiggles will be affected by 
adding a constant power (apart from the too small damp- 
ing coming from b±(k) already discussed), the acoustic 
peak in the two-point correlation function will not expe- 
rience any shift because the one-halo term is so slowly 
varying in Fourier space that its Fourier transform will 
have no structure on 100ft. _1 Mpc scales. In fact, an ex- 
plicit calculation of the acoustic peak shift in the halo 
model was not able to find any shift, placing an upper 
limit of 0.1% [3. 

This argument has been made repeatedly in the lit- 
erature as an explanation of why the acoustic signature 
should be stable under nonlinearities, see e.g. [IS [HI]. 
However, we believe that such arguments are, funda- 
mentally, incorrect. The relevant nonlinear effects are 
not related to virialized objects, but rather to deviations 
from linear perturbation theory at much larger scales (see 
Fig. [lj by the same physics that generates a non-zero bis- 
pectrum at the largest scales (see section IV CI below). 

Furthermore, we now show that the low-fc behavior 
of the one-halo term is in obvious conflict with simula- 



tions. Figure [5] shows the measured mode-coupling power 
Pmc m simulations (symbols with errorbars) as a func- 
tion of scale, for all the large-scale modes available in 
our simulations (corresponding to a fundamental mode 
k = 2tt/1280 /iMpc" 1 « 0.05 /iMpc -1 ). These measure- 
ments were done by measuring the power spectrum and 
the propagator, or cross-correlation between initial and 
final density fields. From Eqs. © and the mode- 
coupling power can be measured in simulations by calcu- 
lating P-G 2 P , that is 



P MC (k) 5 D (k + k') = {5(k)S(k'))- 



fflk)W) 2 
{5o(k)5 (k'))' 



(20) 



where S is the final density field, and So the initial one. 

Note that the mode-coupling power at the largest 
scales is about four orders of magnitude smaller than 
the linear power (e.g. compare to Fig. [I]); we are able to 
measure such small power precisely because of the large 
volume that our simulations probe. The agreement be- 
tween the RPT prediction and simulations is impressive, 
and the failure of the halo model prediction is evident. 
Clearly, using the low-fc limit of the one-halo term to 
make any predictions is, at best, risky. 

We stress that the RPT prediction has no free param- 
eters whatsoever, see Eq. (|2T))) below for an expression 
valid in the low-k limit. As fc — ► 0, the prediction is 



PMc(k,z) 



Pin{q,z) 



d 3 q. 



(21) 



In fact, this low-fc asymptote is a well-known result ex- 
pected from local conservation of momentum: collaps- 
ing structures at small scales can only affect large scales 
through a fc 4 tail [H, H US E3- 

The halo model, on the other hand, predicts that local 
collapse of dark matter particles into halos of size of or- 
der a Mpc can change the dark matter power spectrum 
significantly on Gpc scales! It is not clear to us at present 
how to fix this problem with the one-halo term; although 
one can certainly generate a fc 4 tail by properly including 
mode-coupling physics in the two-halo term, one would, 
in addition, have to somehow cancel the spurious contri- 
bution from the one-halo term. See Section 4.4 in [32j 
for more discussion on this. 

One may wonder whether this discussion is at all rel- 
evant for galaxies, since tracers such as galaxies do not 
conserve momentum because they have significant gravi- 
tational interactions with the dark matter and other trac- 
ers, the fc 4 tail argument does not apply. Thus, a contri- 
bution to the mode-coupling power from virialized halos, 
roughly constant on BAO scales, may be present. Al- 
though given that the low-fc one-halo term prediction is 
in such serious conflict with simulations for the dark mat- 
ter, perhaps one should take specific predictions from it 
for galaxies with a good dose of skepticism. 

Finally we note that the RPT prediction for the mode- 
coupling power is accurate well beyond the low-fc limit. 
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Figure [5] shows the leading contribution to the mode cou- 
pling power in RPT, P M ' C ° P (corresponding to the cou- 
pling of two Fourier modes), and the next order con- 
tribution P m 'q° p (corresponding to the coupling of three 
Fourier modes), which peaks at smaller scales. The sum 
of these two contributions (top solid line) matches the 
simulation results very well up to k ~ 0.2/iMpc" 1 at 
z = 0. For k > 0.2/iMpc _1 three-loop (four-mode cou- 
pling) and higher order contributions not included here 
become important. 



V. SHIFTS 

We now discuss the issue of how nonlinear evolution 
may induce shifts in the acoustic signatures, which can 
potentially bias the determination of the sound horizon 
and lead to systematics in the determination of cosmolog- 
ical parameters such as the equation of state of the dark 
energy. The discussion is based on shifts of the power 
spectrum wiggles and the acoustic peak of the two-point 
function, which one can roughly think of as indicators 
for the size of the sound horizon. A proper assessment 
of the impact of such shifts on the determination of the 
sound horizon depends on the procedure used to analize 
the data, this is beyond the scope of this paper. How- 
ever, our results have something useful to say on what 
such procedures need to take into account. 



A. Power Spectrum Wiggles 

The evolution of the acoustic pattern in the power 
spectrum, Figs. [5] and [3J is rather complex. Looking at 
the position of maxima and minima of the baryon wiggles 
we can notice some difference between linear and nonlin- 
ear cases. While there are definitely shifts between those 
at the percent level, quoting numbers is dependent on 
which quantity one is interested in (see e.g. 14111) . 

Here we follow the method described in [42|], which 
was recently applied to measuring the BAO scale in the 
SDSS and 2dfGRS galaxy power spectra. They define 
a smooth power spectrum out of the observed nonlinear 
power with acoustic oscillations by doing a spline fit to 
a coarse rebinning of it. Then they use the ratio of the 
power to the smooth power and look for oscillation nodes 
in this ratio that oscillates about unity. These nodes are 
then used to calculate the sound horizon scale. 

We define a smooth nonlinear power spectrum by anal- 
ogy with Eq. © 



P s (fc) = G 2 (fc)P s (fc) + P^ c (fc), 



(22) 



do not assume that P MC i s equal to Pmc> i- e - we allow 
for oscillations in Pmc! this is crucial in what follows and 
is why we reach different conclusions than [42] . The ratio 
of nonlinear powers B = P/P s can then be written as, 

B(k) = g(k) B (k) + [1 - g(k)] B MC (k), (23) 

where Bq = Po/Pq is the linear ratio and 

_ G 2 (fc)P s (fc) _ P MC (fc) f0A , 

9( k > = T^m > B MC {k) = 5i 77T - (24) 



P s (fc) 



the idea being that each term is a smooth version of the 
terms in Eq. Note that already the propagator G is a 
smooth quantity, e.g. due to Eqs. (JBHTJ) ; thus Pq is simply 
a smooth version of the linear spectrum Pq. However, we 



To measure the acoustic scale |42j | implicitly assume 
that BMc(k) = 1, in which case it follows from Eq. (123)) 
that the power spectrum nodes defined as -B(fc n odcs) = 1 
coincide with the linear spectrum nodes -Bo(fcnodcs) = 1. 
Thus, they argue, nonlinear evolution does not shift the 
values of fc no d os , which can be related to the acoustic 
scale, e.g. in the approximation of [43| fcnodcs = mr/r s , 
where r s is the comoving sound horizon scale at recom- 
bination. 

However, we now show that the assumption -Bmc(^) = 
1 docs not quite hold, leading to shifts in the values of 
fcnodcs obtained from linear theory. In order to do so, 
it remains to define Pq and P MC appearing in Eq. (|22|) . 
Instead of defining smoothed quantities by rebinning we 
proceed in a slightly different way, for reasons that will 
become clear shortly. The idea, to first approximation, 
is to construct P s out of nonlinear evolution of the linear 
spectrum Pq. 

We define the smooth version of the linear spectrum 
by modifying the BBKS spectrum that approximates the 
overall shape of the wiggly power spectrum with analytic 
functions that are rigorously smooth. A good check is to 
ensure that f PQ(q)dq — f P (q)dq, which in turn will 
lead to a propagator G that is equal for the evolution 
away from smooth and wiggly initial conditions. The 
ratio i?o = Po/Pq so constructed is shown in the left 
panel of Fig. [5] as dashed lines. We then calculate the 
mode coupling contribution for the initial spectrum Pq 
using RPT, let's call this result P^ ooth . This leads to 
a ratio PMc/^Mc°° th that ' s s h° wn m the left panel of 
Fig. [5] as solid lines. 

As it is obvious from this plot, we see that the mode- 
coupling power is not smooth. It contains oscillations 
that are approximately 90 degrees out of phase with the 
acoustic oscillations of the linear power (dashed lines). 
This can be understood from the fact that nonlinear cor- 
rections to the power spectrum are a decreasing func- 
tion of the local spectral index [44]: the dashed curve 
in Fig. [5] modulates the slowly varying spectral index of 
the smooth power, one thus expects maximal (minimal) 
nonlinear corrections at those wavenumbers where the 
derivative of the dashed line is most negative (positive). 
This creates an oscillating correction that is out of phase. 

Also shown in the same plot is the ratio of Pmc mea- 
sured in the simulation (symbols with error bars) ob- 
tained from the measured power and propagator Pmc = 
P — G 2 Pq to the smooth mode-coupling power P^c 00 * 11 
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FIG. 6: Left Panel. Acoustic oscillations are induced in the mode-coupling power spectrum. The solid line denotes the ratio 
of Pmc to Pmc 00 * 11 , with the latter constructed by evolving a smooth spectrum using RPT. The dashed line corresponds to 
the same ratio but in linear theory (Po). The presence of acoustic oscillations in Pmc is also evident in the measurements in 
N-body simulations. Right Panel. The ratios of wiggly to smooth spectra, Po = Po/P) f° r linear theory (dashed lines) and 
Pmc = Pmc/Pmc f° r t ne mode-coupling power in RPT (solid lines) and N-body simulations (symbols with error bars). The 
triangle symbols (with error bars suppressed for clarity) denote Pmc obtained in N-body simulations by the rebinning method. 
The vertical lines denote four linear spectrum oscillation nodes where Po = 1. Note that the oscillations in Pmc are out of 
phase with respect to Po, leading to shifts in the power spectrum nodes toward higher wavenumbers, see Eq. (|25l) . 



calculated from RPT as discussed above. We see that 
the same oscillations are seen in the N-body data, al- 
though there is a significant discrepancy between RPT 
and N-body measurements for k > 0.1 /iMpc -1 , more so 
than in the comparison of total power in Fig. [3J This is 
due to the fact that the measured propagator decays as 
a function of k slightly more slowly than in RPT, thus 
substracting G 2 Po from the measured power leads to a 
lower Pmc- 

To finalize our definition of Pmc m -^q- (j2"2"| we note 

that the ratio Pmc/Pmc^ ( solid line in the left P anel 
of Fig. O does not oscillate about unity as k increases 
due to a small drift. To fix this and thus have a smooth 
nonlinear spectrum about which the wiggly spectrum os- 
cillates, it is simply enough to define P^j C as P^ ooth 
slightly tilted (by 1% at k = 0.2 /i Mpc -1 ) by a linear 
or mild quadratic dependence on k. This brings Pmc in 
Eqs. (|23H24p to a function that oscillates about unity, as 
it should. The same procedure can be done for the N- 
body measurements. The result of this is shown in the 
right panel of Fig. [51 We have thus constructed smooth 
spectra P s about which the wiggly spectra oscillate [54[ . 

We are now ready to see what effect does Pmc (solid 
line in the right panel of Fig. [5J have on the node po- 
sitions determined from P(fc no dcs) = 1 in Eq. (|23p com- 



pared to those in linear theory, Po(fcnodes) = 1 (given 
by the vertical lines in Fig. [5]). Because near the linear 
nodes Pmc is smaller (larger) than unity when Bq has 
positive (negative) derivative, one can easily see graph- 
ically that the net result is to shift all nodes to larger 
wavenumbers. Indeed, close to linear nodes we can write 
B ~ 1 + a (k — feJjJJdS) and Pmc — 1 — e, where the signs 
of a and e always coincide. Then Eq. (|2"3")l says that the 
position of the nodes after nonlinear evolution shifts to 

jU _ j, linear , f ( _^MC_\ . .linear foci 

K-nodcs — "-nodes \ a ) \ Q2ps ) ' ""nodes' l z0 1 

and thus all nodes shift to smaller scales, as expected. 
Note that as k increases, e decreases while a slowly de- 
creases; however, the ratio P^ c /G 2 Pq increases rapidly 
as the mode-coupling power becomes the dominant con- 
tribution to the power (see Fig. [J). Thus, it is difficult a 
priori to estimate how the shift of the nodes depends on 
wavenumbcr. 

Figure [7] shows the results of calculating the 
percentage of shift for the four nodes at k = 
0.057, 0.088, 0.115, 0.14 h Mpc -1 (vertical lines in Fig.© 
as a function of redshift. Note that in this figure we use 
the transfer function for the linear spectrum evaluated at 
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FIG. 7: Percentual shifts due to mode-coupling as a function 
of redshift. The solid line denotes the shift in the acoustic 
peak position of the two-point function. The dashed lines 
denote the shifts in the power spectrum first four nodes in 
the method by f45| (see Fig. |6}. Note that this calculation 
corresponds to dark matter in real space, thus it should be 
considered as a lower limit when redshift space distortions 
and galaxy bias are present. 



redshift z = (unlike the other figures where it was taken 
at z = 49 as for the N-body simulation). The shift as a 
function of redshift for the node corresponding to the fun- 
damental oscillation (node 1) approximately matches the 
shift in the peak of the correlation function due to mode- 
coupling, as discussed in the next section. This should 
not be surprising, the fundamental mode determines the 
location of the peak. 

The shift in the higher-order nodes (2 to 4) depends 
on node number, suggesting that mode-coupling does not 
cause a uniform (independent of k) overall shift. Note 
also that very similar shifts can be derived from our N- 
body results due to the close agreement of Bmg with 
RPT seen in the right panel of Fig. [5] 

These results seem at first glance in disagreement with 
those in || , who using N-body simulations with total vol- 
ume very similar to ours concluded from the dark mat- 
ter real space power spectrum that there were no de- 
tectable shifts. However, many differences in the meth- 
ods used make a direct comparison difficult, e.g. they 
use the rebinning method of [42f to define smooth spec- 
tra, and they do a global fit looking for a uniform shift 
of B with respect to B (going to higher wavenumbers, 
up to k — 0.4 h Mpc -1 ) and simultaneously fitting for 
the damping due to the propagator. A global fit de- 
pends on the behavior of Bmc a lso far from the linear 



nodes, where its effect is smaller, however, is precisely 
there where the details of the damping model, through 
g(k) in Eq. (f2"3"|) , become important. We now try to per- 
form a closer analysis to Their implementation can 
be written as B(k) = \Bn(ka) — l]g(k) + 1, which in view 
of Eq. igSll leads to [Hi, 



B (ka) = B (k) 



1 - g(k) 
g(k) 



[-B M c(*)-l], (26) 



where a is a scaling parameter that characterizes an over- 
all shift in B with respect to Bo. Again, if Bmc were 
equal to unity there would be no overall shift, i.e. a — 1. 
If we use RPT to find the best fit a we obtain, 



1 - a = 0.5%, 0.35%, 0.2% 



(27) 



at z — 0, 0.5, 1 respectively, where we have used scales up 
to k = 0.2/iMpc -1 . 

One source of concern in the method of [42j] is that 
rebinning can introduce noise into the smooth spec- 
trum that could wash out the oscillations in Bmc- We 
have recalculated Bmc from the N-body simulations for 
a smooth power spectrum obtained by rebinning the 
nonlinear wiggly power k 3 P(k) within bins of width 
Afc = 0.044 /iMpc" 1 starting at our fundamental mode 
2tt/1280 w 0.006 /iMpc" 1 . The result is shown as trian- 
gles in Fig. [6j where we have suppressed error bars for 
clarity. We see that rebinning does introduce noise, al- 
though the overall features match what we found using 
our "strictly smooth" method. However, changing the 
binning can affect results. To quantify this, following [f| 
we look for a by minimizing 



(B (k i a)-B (k i )-(?—/pQ) [B MC (h) 
V g(ki) ) 



(28) 

where the sum is over power spectrum bins up to k — 
0.2/iMpc -1 and w t = g{k i )P{k i )/ap{k i ). For the bin- 
ning defined above (Afc = 0.044 /iMpc -1 , triangles in 
Fig. © we find 1 - a = 0.75% ± 0.2% at z = 0, consis- 
tent with our RPT results. If we bin in P(k) instead of 



k i P{k) we get 1- 



0.65%±0.2%. However, trying dif- 



ferent ways of rebinning (e.g. Ak = 0.04 - 0.06 h Mpc" 1 ) 
we found values of 1 — a that range from 1% to 0%, 
which span many standard deviations, suggesting that 
the method of [42| to define smooth spectra is not reli- 
able enough to be used to detect shifts at below the per- 
cent level. These different rebinnings lead to fluctuations 
in the smooth reference spectrum of order 1%, while the 
power spectrum error bars for our total simulation vol- 
ume are of order 0.5% and 0.1% at k = 0.05 h Mpc -1 and 
0.2 /iMpc^ 1 , respectively. 
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FIG. 8: A more detailed view of Figure U close to the acoustic peak. It shows the shift in the position of the acoustic peak 
of the correlation function according to RPT for z = (left panel) and z — 0.5 (right panel). The linear peak (dashed line) 
is smeared out after being convolved with the propagator G, broadening it and shifting its position. The addition of the 
mode-coupling contribution £mc provides another slight shift towards smaller scales, see Fig. [TJ The sum of these two terms 
gives the prediction of RPT (solid line). Vertical dashed lines denote the local maxima of the linear and nonlinear two-point 
correlation function peaks. 



B. Acoustic Peak in the TwoPoint Function 

Figure [5] provides a detailed view close to the acoustic 
peak in the two-point function, showing the linear theory 
(dashed line) and nonlinear (solid) correlation functions, 
together with vertical lines denoting the local peak max- 
ima. The peak shifts towards small scales, and there are 
two sources of shifts. In addition to the smoothing of the 
acoustic peak, convolution by G 2 generates a shift in the 
acoustic peak towards smaller scales, because the linear 
acoustic peak is not quite symmetric about its maximum. 
This shift is referred to as "apparent" in [l|| , and is typ- 
ically taken into account in models used to analyze data, 
e.g. [a. fl3l] . In addition to this, there is another source of 
shifts, termed "physical" in which results from the 
mode-coupling contributions describing coherent infall of 
pairs in the model of [l6| . 

In Figure [7] we separate the contribution to the peak 
shift in these two components, providing the percentual 
shift as a function of redshift due to mode-coupling alone, 
calculated in the one-loop approximation. Note that in 
this figure we have calculated the corresponding shifts 
when the transfer function is taken at z = rather 
than z = 49 (as was used in the N-body simulations). 
This makes the linear peak sharper and thus more robust 
against nonlinearities, decreasing the shift due convolu- 
tion with the propagator by roughly a factor of two, but 



only decreasing the mode-coupling shift by about fifteen 
percent. 

The mode-coupling shift decreases at high redshift, as 
expected. Note that this calculation should be regarded 
as a lower bound to the shifts expected in galaxy redshift 
surveys, which include in addition the effects of redshift 
space distortions and galaxy bias, both of which increase 
their magnitude. See [Ty] for calculations of the acoustic 
peak shift that include the effects of bias [56[ . 

Our results are in broad agreement with the estimates 
of the acoustic peak shift made in [l2| using halof it 
to describe the nonlinear power spectrum. They found 
a shift of 2.4% at z = 0, comparable to ours, e.g. if we 
include the shift due to the propagator at z — we obtain 
1.9%. Our calculations are also in good agreement with 
those in [l6| for z > 0.5, at lower redshifts their method 
overestimates the shifts. As they note, this is due to the 
breakdown of standard perturbation theory. Note that, 
unlike RPT, none of these methods are able to separate 
the total shift into its two different components. 

Our results for the shifts generated by mode-coupling 
are also expected given the results from simulations pre- 
sented in [16| for halos. In that case a fit for the over- 
all shift was attempted in which only the convolution 
("apparent") shift was allowed for. The best fit model 
had statistically significant residuals that made apparent 
that a second source of shifts due to mode-coupling was 



15 



present. Furthermore, note that since this was found for 
point-like halos, the physics responsible for shifts must be 
unrelated to virialized halo profiles. In the next section 
we explain in detail how these shifts arise. 

Finally, note that the shift in the peak of the two-point 
function agrees with that found for the fundamental os- 
cillation node in the power spectrum at all redshifts, as 
expected from the fact that the fundamental oscillation 
sets the scale for the peak of the two-point function. We 
show in the next section that in fact both methods (cor- 
relation function peak and power spectrum oscillation 
nodes) are subject to the same systematics, in the sense 
that if we ignore oscillations in the mode-coupling spec- 
trum the mode-coupling shift in the two-point function 
also vanishes. 



C. Why shifts? 

The shift of the acoustic peak generated by £mc is 
simple to understand from first principles. To do so, 
standard PT (which gives the low-fc behavior of Pmc) 
suffices here [53]. In this case we have 

P MC (k) ~2 J [^ 2 (k-q,q)] 2 P lin (|k-q|)P lin (q)A, 

(29) 

where F 2 denotes the second-order PT kernel (k = k/fc), 

x(Mj-|«u) (30) 

The physical interpretation of these terms is as follows: 
the first term (monopole) describes the growth of pertur- 
bations in the spherical collapse model to second order; 
the second term (dipole) represents the transport of mat- 
ter by the velocity field encoded by the v • V terms in the 
equations of motion; the last term (quadrupole) describes 
the impact of tidal gravitational fields in the growth of 
structure. The excellent agreement between RPT and 
simulations in Fig. [5] at BAO scales tells us that Eq. (j2"9")l 
is quite accurate for the discussion that follows. 

The important contribution for our purposes is the sec- 
ond term in Eq. (|30p . which is the only one that leads 
to derivatives of the linear correlation function. Squar- 
ing the F2 kernel, the lowest order term that leads to a 
shift comes from the double product of the first and sec- 
ond terms, which upon Fourier transformation leads to a 
contribution of the kind, 

&»(r) & } (r), (31) 

where 

tg>(r) =f ■ V"V(r) = J ^PUk)ji(kr) (32) 



is an integrated version of the linear two-point function 
and thus very flat about the acoustic peak. On the other 
hand, the derivative involved in Eq. ([3T|) leads to a net 
shift towards small scales, since £( in is positive (negative) 
to the left (right) of the peak. Close to the peak this 
can be thought as inducing a tilt in the linear correla- 
tion function. Note that shifts appear as a consequence 
of gravitational instability being nonlocal, and by the 
same physics that generates a non-zero bispectrum at the 
largest scales. The spherical model where kernels such as 
F 2 are replaced by their spherical average does not lead 
to shifts of the acoustic peak. Another example that does 
not lead to shifts is the lognormal model [Ty] . 

It's important to note that the two sources of shifts 
discussed above operate in a different way, that due to G 2 
is multiplicative (or convolution) while that due to mode 
coupling is additive. This has interesting consequences. 
The shift generated by G 2 will be absent if the linear peak 
were symmetric, or sharp enough compared to the width 
of G 2 . On the other hand, how much shift mode coupling 
generates depends on how much smoothing G has done: 
since its effect is additive it will be able to modify more 
strongly a weaker peak. Note that these considerations 
apply also in Fourier space. Furthermore, the smoothing 
scale involved in G and the mode-coupling power are not 
very sensitive to features in the linear power spectrum. 
These considerations may help explain why in [9( placing 
a large peak (by BAO standards) in the initial power 
spectrum did not seem to generate shifts. 

Finally, note that Eq. (|2"5)) also explains the main fea- 
tures of Pmc seen in Fig. [S] First, the convolution of 
wiggly spectra leads to oscillations by beating the oscil- 
lation at |k — q| with that at q; note this also leads to a 
result that oscillates not about the mode-coupling power 
of smooth linear spectra, but rather about the average 
of the product of wiggly linear powers, thus explaining 
the drift seen in the solid lines in the left panel of Fig. [5] 
Another way of seeing this, is to calculate Eq. (j2"9")l for a 
power-law power spectrum of arbitrary spectral index n, 
as done in 44]. The result is that Pmc is a strong de- 
creasing function of n: nonlinear growth is enhanced for a 
more negative spectral index. Again, the main cause for 
this effect are the v • V terms that lead to Eq. (|3 1 [) : this 
is also in fact the reason why the skewness of a smoothed 
density field increases as n become more negative. This 
dependence on the local spectral index means that for 
a wiggly linear spectrum, the nonlinear growth will be 
enhanced (suppressed) with respect to the smooth case 
on the nodes of oscillations that have negative (positive) 
derivative, creating the out of phase correction seen in 
Fig. [6j This, in turn, leads to the node shifts derived in 
Eq. (25). 

Since the shift in the peak of the two-point function 
and the power spectrum oscillation nodes have the same 
physical origin, it should not be surprising that switching 
off one of them also makes the other go away. We now 
show that making the mode-coupling power smooth leads 
to no mode-coupling shift for the peak of the correla- 
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FIG. 9: The mode-coupling two-point correlation function as predicted by RPT (solid lines) at redshift z — (left) and z = 1 
(right). The dashed lines show a fit to it of the form given by the second term in Eq. (|34p . where we only allow for an arbitrary 
constant prefactor. The short dash-long dash line denotes the mode-coupling correlation function corresponding to a smooth 
initial power spectrum and therefore has no particular feature at the acoustic scale. Thus, correcting for smooth broadband 
tilts in the power spectrum does not correct for the shift of the two-point function peak. 



tion function. Figure [9] shows (short dashed- long dashed 
line) the Fourier transform of the smooth mode-coupling 
power spectrum obtained by evolving under RPT the 
smooth linear spectrum obtained as discussed in sec- 
tion [V2] We see that the resulting mode-coupling two- 
point function ^ c is flat about the acoustic scale, as ex- 
pected since the smooth initial conditions have no knowl- 
edge of a preferred scale, as opposed to £mc for wiggly 
initial conditions (solid line) which is positive (negative) 
to the left (right) of the linear acoustic scale. For scales 
r < 80 h^ 1 Mpc or r > 130 h^ 1 Mpc, £j^ c approximately 
follows £mCj at levels consistent with the drift seen in 
Bmc about unity in the left panel of Fig. [6j 

From this discussion we see that if the mode-coupling 
power were smooth (5mc = 1), there would no shift 
in the peak of the two-point function caused by mode- 
coupling. In other words, correcting for smooth broad- 
band tilts in the power spectrum does not correct for the 
shift of the two-point function peak. 



D. A Simple Model 

We can now estimate how the shifts generated by mode 
coupling depend on galaxy bias and redshift distortions. 
Indeed, in the presence of local bias perturbation theory 



tells us that Eq. gets modified into [3 HE [H[ 

P^c(k) = 2M 2 J F 2 (k-< l ,<i)P ii n(\k-c i \)P lin (q)d 3 q 

+ b\P uc {k) + ^ J P lin (|k-q|)P lin (q)d 3 9 , 

(33) 

where b\ and b 2 are the linear and quadratic bias param- 
eters. We see that an additional term that induces shift 
in the acoustic peak of the two point function arises from 
nonlinear bias, proportional to b\b 2 . It has the same form 
as Eq. (|31[) and will lead to enhanced shifts compared to 
the linear bias case for galaxies that preferentially popu- 
late high-mass halos, in which case b 2 > 0, in agreement 
with the arguments given in [14j for the power spectrum. 
In addition, redshift-space distortions amplify the effect 
of the second term in the F 2 kernel responsible for shifts 
by (1 + /3/i 2 ) [H, where (5 = f/h, f = din D + / din a, 
and \i is the cosine of the wave- vector along the line of 
sight, under the approximation that the mapping from 
real to redshift space can be treated perturbatively. 

Motivated by the above discussion we can suggest a 
simple phcnomcnological model to take into account non- 
linear effects close to the acoustic peak that may be useful 
in practical applications. One can try to fit the observed 
galaxy correlation function monopole in redshift space by 
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Us(r) = A [e- r2 ^ 2 ® £ lin ](r) +B^(t) (34) 

where we have used that ^ in Eq. flUT) varies very lit- 
tle around the acoustic peak, and A, B and a are con- 
stants that can be fit for. From linear bias and red- 
shift distortions, one expects A ~ b\ (1 + 2/3/3 + /3 2 /5), 

B ~ ^linOv) < A and a 2 - 4<t 2 (1 + //3), where r p 
is the position of the acoustic peak. Similarly, in the 
method that uses the power spectrum oscillation nodes, 
one may want to add a small correction that has an out of 
phase oscillation, which is analogous to the second term 
in Eq. (|34[) . Finally, note that these correction terms 
coming from mode coupling are intrinsically second or- 
der in the power spectrum or correlation function, thus 
one may want to include this information in the model- 
ing. 

It is beyond the scope of this paper to study such mod- 
els in detail, but for the same reason that a result from 
first principles such as Eq. ^ was later shown to be use- 
ful phenomenologically in the presence of bias and red- 
shift distortions (see e.g. @,[l3j]) we would like to stress 
that the second term in Eq. (|34|) . despite being derived 
here invoking a number of approximations, is likely to 
prove similarly useful. To see this, we show in Fig. [5] the 
full result for £mc in RPT and a fit to it of the form given 
by the second term in Eq. (|34[) . where we only allow for 
an arbitrary constant prefactor. We see that despite the 
simplicity of the model it captures the behavior rather 
well. One could certainly improve the fit by including 
more terms that result from Eq. (J21J). 

In addition, the relevance of Eq. (j3~4"|) for more compli- 
cated situations where bias and redshift distortions are 
present can be appreciated from looking at the form of 
the residuals to fits of halo-halo correlation functions in 
redshift space (Fig. 5 in [l6j]) that includes only the first 
term in Eq. (|34|) . These residuals have a form very well 
matched by the second term in Eq. (f3~4"|) , shown in Fig. [5] 
by dashed lines. 



VI. CONCLUSIONS 

We calculated the nonlinear power spectrum and two- 
point correlation function in renormalized perturbation 
theory ([H, HI], RPT), with particular emphasis on the 
description of baryon acoustic oscillations (BAO). Our 
predictions are in excellent agreement with detailed mea- 
surements in numerical simulations, as shown in Figs. [5] 
through [8] 

In this paper we have concentrated on a single cos- 
mological model at different redshifts; this was driven 
by the desire to compare our predictions to simula- 
tions at the percent level at BAO scales, which requires 
very large volumes, e.g. our total simulation volume is 
~ 105(/i _1 Gpc) 3 . Calculating RPT predictions for dif- 
ferent cosmological models is rather simple, while having 



corresponding simulations with the required precision at 
BAO scales is significantly time consuming, this will be 
considered in a forthcoming paper (23[. 

RPT describes the power spectrum (or correlation 
function) as the sum of two contributions, Eq. ©, both 
of which are sensitive to nonlinear effects. One term dom- 
inates at large scales and is proportional (via the propa- 
gator G) to the initial power spectrum at the same scale, 
reducing to linear theory in the large-scale limit. Wc 
already tested in [lj| the RPT prediction for the prop- 
agator against numerical simulations, for both densities 
and velocities as a function of redshift. The propagator 
is responsible for the damping of BAO in the power spec- 
trum, or smoothing of the acoustic peak in the two-point 
function. 

The second contribution, the mode-coupling power 
spectrum Pmc, describes the generation of power at a 
given scale due to the coupling of two or more modes from 
different scales. Calculation of this in RPT was the main 
focus of this paper (see Appendix [Alfor technical details). 
Figure [5] shows that the mode-coupling power in the sim- 
ulations grows as dictated by RPT. Mode coupling modi- 
fies the acoustic pattern in the power spectrum, inducing 
oscillations out of phase with those in linear theory, and 
this generates a shift towards small scales in the acoustic 
peak of the two-point function. Our predictions for these 
effects are also in excellent agreement with simulations. 

The halo model describes two-point statistics in a way 
that resembles RPT. However, there are important dif- 
ferences that we highlighted in this paper. Although the 
sum of the one-halo and two-halo terms as written in 
section IIVI describes the power spectrum in simulations 
to within ~ 15%, this is highly misleading because the 
two-halo term is too large and the one-halo term too 
small in the range of scales between linear and virialized 
regimes, 0.05 /iMpc" 1 < k < l/iMpc^ 1 . This is simply 
because in this range the physics that operates is nei- 
ther linear theory nor due to virialized halos. This is 
obvious when one tries to use the halo model to describe 
higher-order statistics such as the bispectrum: using lin- 
ear theory plus virialized halos fails miserably, and one 
must include weakly nonlinear physics [4?| . Once this 
is recognized, it is simply not optional to include or not 
such contributions for two-point statistics. Similarly, as 
emphasized in [14], one must pay attention to nonlin- 
earities in the large-scale bias relation when describing 
two-point statistics of halos or galaxies at large scales. 

Fixing the two-halo term by adding mode-coupling 
contributions and correct damping of the linear spectrum 
is possible by incorporating the missing physics (e.g. by 
RPT); however, it is less clear at present how to avoid the 
wrong low-fc limit of the one- halo term (see Fig. [S]). In 
any case, the test any improved model must pass is very 
simple: it must match the cross-correlation of initial and 
final density fields and the mode-coupling power, instead 
of only using the power spectrum or correlation function 
as a diagnostic. 

Having demonstrated that the mode-coupling power 
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predicted by RPT is in very good agreement with sim- 
ulations, we discussed what are the implications for the 
evolution of the BAO. We showed that the mode-coupling 
power contains oscillations that are essentially out of 
phase with those in the linear spectrum. This leads to 
shifts towards smaller scales in the position of the os- 
cillation nodes used in the BAO analysis method pro- 
posed in [42j. We showed that this mode-coupling con- 
tribution, when Fourier transformed, naturally leads to 
a shift of the acoustic peak in the two-point function to- 
wards smaller scales, through terms proportional to the 
derivative of the linear correlation function, Eq. (f3"Tj) . We 
also showed that correcting for smooth broadband tilts 
in the power spectrum does not correct for the shift of 
the two-point function peak. 

We gave predictions for how these mode-coupling shifts 
depend on redshift, which should be considered as a ro- 
bust lower limit to the more realistic case when galaxy 
bias and redshift distortions are also present. See [lj| for 
calculations that also include the effects of galaxy bias 
on the acoustic peak shift. 

Once the source of shifts is identified, however, it is 
easy to see what to expect when redshift-space distor- 
tions and galaxy bias are included. This understanding 
allows us to suggest a simple physically motivated model 
that can be used in dealing with such effects in observa- 
tions. In fact, the model naturally predicts the kind of 
pattern observed in ,16| as residuals when the halo cor- 
relation function in redshift space is fit to a model that 
only includes the effects of smoothing of the peak. 

Finally, we would like to emphasize that the nonlincar- 
ities we focus on in this paper do operate at the acoustic 
scale. The same mode-coupling effects that we discussed 
here are responsible for generating a non-zero bispectrum 
at large scales; e.g. this is the reason why the skewness in 
spheres of diameter equal to the size of the acoustic scale 
is about forty percent at z — [48]. Most importantly, 
the detailed form of these couplings, e.g. Eq. (j30|) . has 
been observed in galaxy surveys, giving rise to the depen- 
dence of the bispectrum on triangle shape (49j . It might 
be possible, in fact, to empirically determine corrections 
for clustering systematics in two-point statistics from the 
measured three-point function. 
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APPENDIX A: CALCULATING THE POWER 
SPECTRUM IN RPT 

The nonlinear evolution of the power spectrum in the 
framework of RPT is made out of two basic ingredients: 
the nonlinear propagator G and the mode coupling power 
Pmc UM- Recalling Eq. @ we have for the density field 
power spectrum, 



P(k, z) = G 2 (k, z) P Q {k) + P MC {k, z), 



(Al) 



where Po(k) is the initial spectrum of fluctuations. In 
[l|| we computed the propagator G while in this paper 
we completed the description by computing Pmc at the 
scales relevant for BAO. Although a detailed analysis of 
these calculations will be given elsewhere [23| we present 
in this appendix the expression for the contributions to 
Pmc that we use throughout this paper (e.g. Eq. (flQ|) ). 
together with the approximations involved. 

In [l8| we showed that Pmc can be written as an infi- 
nite sum of contributions expressed in terms of only three 
resummed quantities, the nonlinear propagator G a b{k, z), 
the full vertex T a b c (k, ki, k2) and the nonlinear power 
spectrum P a b(k,z), where sub-indices here each take on 
two values since they refer to either density contrast or 
velocity divergence fields. To lowest order these quanti- 
ties reduce to the linear propagator (77 = lnD + ), 



9ab(z) 



3 2 
3 2 



-3rj/2 



-2 2 
3 -3 



(A2) 



which one can think of as four growth factors, i.e. how 
much does the density or velocity grow with respect to 
the initial densities and velocities (thus making a two 
by two matrix) and, similarly, four decaying modes; the 
symmetrized vertex, 
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with ki2 = ki+k 2 and ^ a bc — otherwise; and the linear 
power spectra 



P^(k,z) = [D + (z)¥P (k) 



1 1 
1 1 



(A6) 



from standard PT [l8|, [2J, |50|. Note that we are as- 
suming growing mode initial conditions where the initial 
density contrast S and velocity divergence V • u (where 
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v = —fHu, with / = d In D + /d In a, H = dlna/dr and 
t conformal time) are the same random field. This is the 
reason why all the linear spectra in Eq. (|A6j) are deter- 
mined by a single spectrum, Pq. For the same reason, 
the combination G a = G a b ■ (1,1) represents the prop- 
agator for density and velocities for growing mode ini- 
tial conditions, and it appears frequently. For instance, 
the density propagator in Eq. (|Aip is actually given by 
G = G\ = Gn + G12. 

As we mentioned in Sec. IIII Bl each contribution to Pmc 
is dominant only in a narrow range of scales centered at 
an increasingly higher value of fc. We found that consid- 
ering the first two contributions to Pmc was sufficient to 
describe the range of scales relevant to BAO (see Figs. [3] 
and EI). Since in diagrammatic language they correspond 
to a one-loop diagram and a two-loop diagram [18| we 
refer to them as P M c° P and P M c° P , respectively. 

In (l9l | we studied in detail the nonlinear propagator 
G a b(k,z) and derived an analytical prescription for it, 
valid for all times and scales, that showed a remarkably 
good agreement with measurements in numerical simu- 
lations well into the nonlinear regime (see Sec. IIII Bl for 
a brief account of the main ideas and procedure behind 
this derivation). Using better simulations, we have found 
that our prediction was slightly overestimating the new 
measurements by a few percent at intermediate scales. 
Out of the few reasons that could cause this slight mis- 
match the leading possibility is a small contribution from 
subdominant diagrams in the large-fc limit resummation 
for G a b- Research under way will confirm whether this 
is the reason or not, for this paper we proceed to make 
a small correction to our calculated propagator in [l9j 
under this assumption. It is important to note that once 
this is done, the prediction for the nonlinear power spec- 
trum is fully determined. As we discussed in connection 
with Fig. [51 the simulations used in this work have a 
propagator that decays slightly slower than our modified 



RPT prescription; that is, our procedure is not a fit to 
the propagator measured in the simulations. 

As mentioned in Sec. IIII Bl the prediction for G a b was 
obtained by matching its low-fc behavior, computable ex- 
actly using one-loop PT, with its large-fc limit, where it 
can be shown to decay exponentially as exp(— fc 2 tr^/2) 
with al = |J P L (k,z)d 3 q/q 2 (see Eq. ©). In turn, this 
limit was obtained by resumming an infinite sub-set of 
diagrams in the series expansion for the propagator. In 
.18]" we argued that this sub-set gives the dominant con- 
tribution to the resummed propagator. Although this 
remains true we have estimated the contribution com- 
ing from sub-leading diagrams in the resummation of the 
large-fc limit and found that it can reproduce the few per- 
cent level disagreement with the simulations at all red- 
shifts. 

The way we estimated this subdominant contribution 
is as follows. We notice that a sub-set of these sublead- 
ing diagrams leads to power spectrum resummation. In 
this case the nonlinear spectrum, instead of the linear, 
should be used to compute the nonlinear propagator in 
its large-fc limit (i.e. a 2 — > a 2 {z)a 2 , with a 2 (z) given by 
/ P n \(k,z)d 3 q/q 2 / J Pn n (k, z)d 3 q/q 2 ). For simplicity, we 
computed this factor using halof it to describe P n i(fc, z), 
and found that it grows monotically from 1 at high red- 
shift to a ~ 1.05 at z = 0. Thus we multiply by a 2 (z) the 
exponents in each component of G a b in Eq. (41) of [l9T |. 

For the present work, we use the functional form for 
Gab given in Eqs. (41) of [19], corrected as described 
above, to compute P^^ and Pmc° P - ^ n addition, the 
full vertex is set at its lowest order, i.e. the symmetrized 
vertex is given by Eqs. (|A3flA5|) . while the nonlinear spec- 
trum that enters into the calculation of Pmc is given by 
P a b(k,z) = G a (k, z)P (k)Gb(k 1 z), which in a iterative 
scheme based on Eq. (I All) would correspond to the first 
step. Under this approximation the one-loop contribu- 
tion to Pmc is given by (see Eq. (33) in [Hj]), 



with 



PiLoop(k,??) = 2 J d 3 qAiLoo P (k,q,77) AiLoo P (-k,-q,77)Po(g)P (|k-q|), (A7) 

A 1Loop (k,q,77) = / d Sl G 16 (|k|,7 y , Sl )7^(k,q,k-q)G c (q, S i,0)G rf (|k-q|, Sl ,0), (A8) 
Jo 

where the time variable rj — lnP + , and G a = G a \ + G a2 . The second contribution, corresponding to a two- loop 
diagram, is given by (see Eq. (34) in [ill]). 

P 2 Loo P (k,7 ? ) = wjd 3 q yd 3 pA 2Loop (k,q,p,77) A 2Loop (-k,-p,-q,77)Po( g )Po(p)Po(|k-p-q|), (A9) 



with, 



A 2L oo P (k,q,p,77) = f dsi f ds 2 Gi h (|k|, 77, s x ) j^iK k - q, q) G c/ (|k - q|, s x , s 2 ) G d (q, s 1; 0) 
Jo Jo 

7 g h (k - q, p, k p q) G g (\p\, 82, 0) G fc (|k p q|, S2, 0). (A10) 
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These lowest order assumptions for the full vertex and 
nonlinear power in the computation of Pmc are expected 
to be very robust at low-fc, e.g. if one were to stop at the 
one-loop level where the dominant behavior is dictated 
by the decay of the propagator. At the two-loop level 
one should also check the impact on P^l'^' p coming from 
the renormalization of the vertex or a second step in the 
iterative scheme. At present we are studying these effects 
and the results will be presented in [23| and . 

Nonetheless it can be concluded a posteriori that this 
lowest order approach for T and P works very well at 
the BAO scales given the excellent agreement shown in 
Figs. rjKSJ However, we expect that at the one-percent 
level, some of our approximations within the framework 
of RPT are likely to begin to fail. We can briefly sum- 
marize the main issues at play, 

i) The matching ansatz for the propagator done 
in and the size of sub-leading diagrams in the 
large- k limit. 

ii) The so-called "vertex renormalization" [l8| and the 
inclusion of mode coupling power in the r.h.s of 
Eqs. JATlElO) mu st be quantified. 

iii) The incorrect cosmological dependence of decay- 
ing modes is expected to play a role in the nonlin- 
ear regime, as perturbations do not remain in the 
growing mode after many interactions. This enters 
through our assumption that the linear decaying 
mode obeys D_ oc D+ 3 ^ 2 , see Eq. (|A2|) . This is 

only true as long as we can set //f2„ — to 
unity 

iv) The impact of shell crossing, which so far we have 
neglected since we work in the fluid approximation. 

The first set of issues were effectively corrected by 
renormalizing a v as described above. But further re- 
search is needed to incorporate this, if needed, in a 



stricter way. The three remaining effects should have 
negligible impact on the modeling of BAO since they are 
expected to become important at smaller scales, but they 
are approachable and will be discussed in future work. 

Regarding the last issue listed above, it was recently 
argued in [513 that the impact of shell-crossing on the 
power spectrum at BAO scales is of order one percent, 
and that any agreement at the percent level (such as 
that shown in Fig. [3|) must be "accidental" unless justi- 
fied by solving the collisionlcss Boltzmann equation. In 
principle percent-level corrections from shell crossing are 
possible, although we disagree that the calculation in [5l| 
addresses the issue. The impact of shell crossing must be 
estimated from the difference between the collisionless so- 
lution and the fluid approximation. In [5lT ] , on the other 
hand, the argument is based on a comparison between 
two extremely different situations (and none of them is 
the fluid limit): the collisionless case versus the "sticky" 
dark matter scenario where collisions have a huge cross 
section. It is not clear that the difference between these 
two situations would be a good estimate of the effect of 
shell crossing (i.e. the fluid limit may be a much bet- 
ter approximation in the collisionlcss case). Most im- 
portantly, the halo model as written in section IIVI was 
used in place of the collisionless solution, and the halo 
model with delta function profiles (but same mass func- 
tion) was assumed to describe the collisional solution. As 
discussed in section ITVl the one-halo term is about a fac- 
tor of four too small at k ~ 0.1 /iMpc -1 (see Fig. [5j, and 
consequently the two-halo term (where most of the effect 
quoted in [5 11 ] is coming from) is too large by about the 
same factor. Given all these assumptions, we do not feel 
the conclusions in [5l| are warranted. 

One can estimate directly from Fig. [5] when the con- 
tribution from virialized profiles (where shell crossing 
has definitely taken place) becomes comparable to the 
mode-coupling power in our simulations by extrapolat- 
ing the one-halo term to small scales, and this happens 
at k ~ 1/iMpc" 1 . 
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